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Abstract 

We study an inhomogeneous sandpile model in which two different toppling rules are defined. 
For any site only one rule is applied corresponding to either the Bak, Tang and Wiesenfeld model 
[P.Bak, C. Tang, and K. Wiesenfeld, Phys. Rev. Lett. 59, 381 (1987)] or the Manna two-state 
sandpile model [S. S. Manna, J. Phys. A 24, L363 (1991)]. A parameter c is introduced which 
describes a density of sites which are randomly deployed and where the stochastic Manna rules are 
applied. The results show that the avalanche area exponent r a , avalanche size exponent r s , and 
capacity fractal dimension D s depend on the density c. A crossover from multifractal scaling of 
the Bak, Tang, and Wiesenfeld model (c = 0) to finite size scaling was found. The critical density 
c is found to be in the interval < c < 0.01. These results demonstrate that local dynamical rules 
are important and can change the global properties of the model. 
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I. INTRODUCTION 

Bak, Tang, and Wiesenfeld (BTW) [1] introduced a concept of self-organized criticality 
(SOC) as a common feature of different dynamical systems where the power-law temporal 
or spatial correlations are extended over several decades. Dynamical systems with many 
interacting degrees of freedom and with short range couplings naturally evolve into a critical 
state through a self-organized process. They proposed a simple cellular automaton with de- 
terministic rules, which is known as a sandpile model, to demonstrate this new phenomenon. 
In this model the relaxation rules are conservative, no dissipation takes place during relax- 

n 

ation, and correspond to a nonlinear diffusion equation [1]. Generally, the sandpile model 
is represented by a <i-dimensional hypercube of the finite linear size L. Its boundaries are 
open and allow an energy dissipation, which takes place only at the boundaries. 

Manna proposed a two-state version of the sandpile model where no more than one 
particle is allowed to be at a site in the stationary state. If one particle is added to a 
randomly chosen site, then relaxation starts depending on the occupancy of the site. If the 
site is empty, a particle is launched. In the case when the site is not empty, a hard core 
interaction throws the particles out from the site and the particles are redistributed in a 
random manner among its neighbours. All sites affected by this redistribution create an 
avalanche. An avalanche is stopped if any site reached the stationary state, i.e. no more 
than one particle occupies a site. 

The first systematic study of scaling properties, universality and classification of determin- 
istic sandpile models was carried out by Kadanoff et al. [3j]. Using numerical simulations and 
by varying the underlying microscopic rules which describe how an avalanche is generated 
they investigated whether different models have the same universal properties. Applying 
finite-size scaling (FSS) and multifractal scaling techniques they studied how a finite-size of 
the system affects scaling properties. 

The real-space renormalization group calculations Q] suggested that deterministic [l| 
and stochastic sandpile models belong to the same universality class. On the other hand, 
many numerical results ^, O, Q, [si, [^j show clearly two different universality classes. They 
do not confirm the hypothesis that small modifications in the dynamical rules of the models 
do not change the universality class, presented by Chessa at al. [lo| . 



This study was motivated by the results published by Tebaldi et al. [llf, and Stella 
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and Menech [12j], where a multifractal scaling of an avalanche size distribution of the BTW 
model was demonstrated. They assume that a multifractal character for SOC models like 
the BTW model is a crucial step towards the solution of universality issues. By applying 



the moment analysis they found FSS for the two-state Manna model 12J. Based on these 



results they conclude that the 2D BTW model and the Manna model belong to qualitatively 



different universality classes. This assumption was confirmed recently [13l . Il4j |. where a 
precise toppling balance has been investigated in more detail. 

In this paper we report the results of disturbing the dynamics of the BTW model using 
stochastic Manna sites which are randomly deployed. They can introduce stochastic events 
during an avalanche propagation. Our model was derived from the inhomogeneous sandpile 



model [15|] in witch two different deterministic toppling rules were defined. In the proposed 
model the first toppling rule corresponds to the BTW model [l| and the second rule is now 
stochastic and corresponds to the two-state Manna model [21. The model is similar to that 

n 

in Ref. |14j], however we applied the original toppling rules of the listed sandpile models. 

The paper is organized as follows. The inhomogeneous sandpile model is introduced in 
Sec. [Ill The avalanche scaling exponents, capacity fractal dimensions and crossover from 
multifractal to FSS are investigated with numerical simulations and the results are presented 
in Sec. IIHI The Sec. [IV] is devoted to a discussion which is followed by conclusions in Sec. 

El 

II. MATHEMATICAL MODEL 

We consider a d- dimensional hyper-cubic lattice of linear size L, and a notation presented 
by Ben-Hur et al. [7] is followed to define a sandpile model. Each site i has assigned a 
dynamical variable E(i) that generally represents a physical quantity such as energy, grain 
density, stress, etc. A configuration {E(i)} is classified as stable if for all sites E(i) < E c , 



where E c is a threshold value. We note that the two-state Manna model [2| has no threshold 
E c . The Manna model has defined a hard core repulsion interaction among different particles 
at the same position. This hard core repulsion interaction can be described by a threshold 
where the threshold value E c = 2 is assigned to any site. In our inhomogeneous sandpile 
model, the threshold values E c depend on the site position i, E c (i) jTsj] . The conditions for a 
stationary state, a stable configuration {E(i)} (no avalanche), are now E(i) < E c (i), where 
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the threshold E c (i) at the site i was randomly chosen from two allowed values 

E c (i)={ (1) 
[e» = 2. 

For any site i the threshold E c (i) [Eq. ([T])] is defined in such a manner that n randomly 
chosen sites have the value E]. 1 and the remaining L d — n sites have the value E c . The 
density of sites with the threshold value E]} is denoted c, and c = n/L d . 

Let us assume that a stable configuration {i?(j)} is given, and then we select a site i 
at random and increase E(i) by some amount 5E. We now consider 5E = 1 for any site. 
When an unstable configuration is reached, E(i) > E c (i), a relaxation takes place. An 
unstable site i lowers its energy, that is distributed among the neighbor sites. The directions 
to the neighbor sites are defined by the vectors ei = (0, 1), e 2 = (0, —1), e 3 = (1,0), and 
e 4 = (—1,0). The relaxation is defined by the following rules 

£(i)^£(i)-5>£(i), (2) 

e 

E(i + e) ^£(i + e) + A£(e), (3) 
^AE(e) = E e (i), (4) 



e 



(5) 



{ei, e 2 , e 3 , e 4 } if E c (i) = E{ 
{e c , e v } if E c (i) = Ei 1 

where e is a set of vectors from the site i to its neighbors. The indexes ( and rj are integers 
1, 2, 3, and 4 randomly chosen at any relaxation. The neighbors that receive the energy can 
became unstable and topple, thus generating an avalanche. The distribution of energy is 
described by Eqs. and ©, we added additional rules Eqs. (jlj) and © which specify 
the manner how the energy is distributed depending on the position i, threshold E c {\) [Eq. 
PJ], and corresponding sandpile model [Eq. (JHJ)]. The relaxation rules Eqs. (J2J)- (EJ) are 
applied until that moment when a new stable configuration is reached again, for all sites 
^(i) < E c (i). Obviously, during one avalanche an arbitrary unstable site i can transfer the 
energy E c (i) a few times to became stable, E(i) < E c (i). A <i-dimensional lattice has open 



boundaries so added energy can flow outside the system, and an energy dissipation takes 
place only at the boundaries. 

This model has been designed to enable a well defined change between two well known 
nondirected sandpile models: deterministic [l| and stochastic j^J (nondirected only on av- 
erage) similarly as in Ref. [3]. The model belongs to the critical height models with 
conservative relaxation rules and with undirected energy transfer where the two thresholds 
are randomly frozen. It can be characterized as a sandpile with a possibility to modify its 
scaling behaviors. 

III. RESULTS 

We shall report the results obtained using numerical simulation of the conservative, undi- 
rected, critical height sandpile model defined by Eqs. (EJ)- (J5I) - The simulations were carried 
out for the following parameters: d = 2, two-dimensional lattice of linear sizes L = 256, 512 
and 1024, randomly added energy 5E = 1, two thresholds either E T C = 4 or E* 1 = 2, and 
with density of sites with threshold E], 1 in the interval < c < 1. In our simulations we 
have used the density c as a model parameter. For densities of stochastic sites c = and 
1 the model behaves as the BTW model [l| and Manna model [2J, respectively, which are 
both considered to be Abelian [1 61 ] . 

Avalanches can be characterized by such properties as their size, area, lifetime, linear 
size, and perimeter. We concentrate only on a minimal number of parameters which are 
necessary to demonstrate the investigated phenomena: the avalanche area a and avalanche 
size s. Here the avalanche area a is the number of lattice sites that have relaxed at least once 
during the avalanche. The avalanche size s is the total number of relaxations that occurred 
during the avalanche. The probability distributions of these variables are usually described 
as power-laws with cutoff 

P{x) = x- Tx E{x/x c ), (6) 

where x = a,s. When the system size L goes to infinity, the cutoff x c diverges as x c ~ L Dx . 
If we assume FSS, then the set of exponents (t x ,D x ) from Eq. ([6]) defines the universality 
class of the model [lol ]. 

The avalanche area probability distribution P(a) and avalanche size probability distribu- 
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tions P(s) have been analyzed at finite lattice sizes L = 256,512, and 1024. It is expected 
that these distributions follow a power-law P(x) ~ x Tx [Eq. (jSJ)]. For any lattice size L 
and density c the corresponding scaling exponents t x ,l(c) were determined. The scaling 
exponents found in the numerical simulations for the largest lattice size L = 1024 and for 
selected densities c are presented in Table [H It is evident that the exponents are increasing 
with c in the interval < c < 0.1 and then for densities c > 0.1 they are almost constant. 

The scaling exponents t Xj l show a finite size-effect when the lattice size L is changed. 
Their dependences on lattice sizes L are approximated by a formula proposed by Manna 



const. 

x = x L ^ - (7) 

This approximation was used to extrapolate the scaling exponents T Xt i^ fOC for the infinite 
lattice L —>■ oo. 

The avalanche size probability distributions P(s) obey the power-law dependence for any 
density c. The corresponding scaling exponents t Si l-^oo(c) are shown in the Fig. CD In the 
range of densities 0.01 < c < 0.1 these scaling exponents decrease from ^^^^(O.Ol) = 
1.37±0.025 to ^^^^(O.l) = 1.29 ± 0.025 and then, for higher densities c > 0.1, are almost 
constant. 

The avalanche area scaling exponents T a L^oo show a more complex dependence on the 
density c. For densities 0.09 < c < 0.5 they decrease from 7^^00(0.09) = 1.49±0.025 to 
T a ,L^oo(0.5) = 1.38 ± 0.025, then for higher densities c > 0.5 the exponents r aj i^oo(c) are 
almost constant. It was observed that for densities 0.01 < c < 0.09 the avalanche area 
distributions P{a) do not follow exactly a power-law dependence as it is expected from Eq. 
(E]). Therefore the exponents ^ ^^^(c) from this density interval are not included in Fig. 
[TJ One typical example is shown in Fig. [5] where the density of random toppling sites is 



TABLE I: The scaling exponents t X) l=io2a{c) for the finite lattice size L = 1024 and selected 
densities in the interval < c < 1. The statistical errors are ±0.001. 

density c 0.01 0.10 0.50 1 

r a ,L=i024 1-131 - 1.291 1.315 1.338 

t s ,l=io24 1-137 1.240 1.263 1.266 1.283 
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density c 



FIG. 1: (Color online) The avalanche area and size scaling exponents T a L-*oo an d t S) l-^oo are 
approximated for the infinite lattice size L — > oo. The exponents depend on the density c of 
Manna sites. 




avalanche area a 



FIG. 2: (Color online) The avalanche area distribution P(a) does not follow exactly a power-law 
function. The parameters used in the numerical simulation were: density c = 0.01 and linear lattice 
size L = 1024. 

c = 0.01 and the lattice size is L = 1024. The double-log plot of area distribution function 
P{a) clearly shows that a possible approximation function is not a straight line which must 
correspond to the simple power-law dependence. 

For the two well known sandpile models, BTW (c = 0) and Manna (c = 1) the scaling 
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FIG. 3: (Color online) The extrapolated spectra of the avalanche area a and avalanche size s for 
various densities of Manna toppling sites: (a) c = BTW model, which shows multifractal scaling, 
(b) c = 0.01 at which the multifractal scaling of BTW model is destroyed, and (c) c = 0.95 where 
the model shows the FSS near the two-state Manna model (c = 1). The maximal error bars of 
f(a,c) are for q ~ 0, and are approximately ±0.05, but for a higher q they are smaller. The a 
values are determined within errors ±0.025. 

exponents ^^(O) = 1.26, ^^(O) = 1.23, ^^(l) = 1.36, and ^^(l) = 1.27 
were found. In addition, for all densities c (see Fig. CP the relation r aj L-*oo{c) > Ts,l^oo(c) 
is valid. 

The scaling exponents t Xj l as functions of the lattice size L show a finite-size scaling effect 
[Eq. ((7j)]. An exact determination of scaling exponents r^^oo from numerical experiments 



S 



is therefore a difficult task. A new method was introduced [6| to increase the numerical 
accuracy of the exponents based on their direct determination. We found that the method 
gives slightly larger exponents than a simple extrapolation of Eq. (J7J). However, the expo- 
nents t s do not fluctuate around their mean values as it was observed in the paper p. Our 
error bars were larger, therefore we have to repeat this analysis again in more details. 



Tebaldi et al. 



ll| found that in the BTW model the avalanche area distributions P{a) 



show FSS and avalanche size distribution P(s) scale as a multifractal. To describe these 
scaling properties rather a multifractal spectrum f(a) versus a than the single scaling ex- 
ponent t s [Eq. (jHJ)] is necessary. Thus, the scaling exponent r s loses the importance and is 
replaced by a spectrum of exponents. Despite this fact, the avalanche size scaling exponents 
Ts,l^oo(0) are determined. They enable a comparison with the previous results, since the 
whole point is that the exponent ^ ^^^(O) does not exist. The recent studies [ll|, ll2| led us 
to analyze the multifractal properties of the model given by Eqs. ©-© for various densities 
c. To determine the multifractal spectra a method presented in the paper [12j was useful. 
There, for any finite-size lattices L, the quantities a x (q,L) = (log(x)x q ) / [log(L) (x q )} and 
a x (q, L) ~ log {{x q )) /log(L) were computed. It was observed that a x (q, L) and o~ x (q, L) show 
a finite-size dependence on the system size L, which is well approximated by Eq. (j7]) and 
this relation was used to extrapolate L — > oo quantities. Based on the Legendre structure 
relating f x to a x , a parametric representation of f x (a x ) by plotting f x (q) = cr x (q) — a x (q)q 
versus a x (q) can be obtained [l^ . 

Some significant spectra of f x (a x ,c) extrapolated for an infinite lattice size L — > oo are 
shown for illustration in Fig. [31 The f x (a x ,c) values were determined for the parameter 
q in the range —3.5 < q < 3.5 and they are limited by errors about ±0.08, similarly as in 



Ref. 



121 ] . We have observed that if f x (a x , c) spectra are computed for all avalanches where 



a > 50 then the errors of f x (a x , c) are ±0.05. The multifractal scaling of the avalanche size 
probability distribution P(s) and FSS of avalanche area probability distribution P(a) were 
found at density c = (see Fig. [3] (a)). The avalanche probability distributions P(x) show 
FSS for densities c = 0.01 Fig. 0(b), and for c = 0.95 Fig. [3](c) which is close to the Manna 
model (c = 1). The spectra for c = and 1 agree well with the previous results [12]. It was 
found that the multifractal scaling of P(s) was destroyed (Fig. E^b)) at a relatively small 
density of Manna sites < c < 0.01. 

Stella et al. 12] claim that if probability distributions P(x) satisfied FSS the large q data 
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FIG. 4: (Color online) A crossover from multifractal scaling to finite size scaling takes place at 
< c < 0.01. The hatched area border an interval of A/(c) in which A/(c) = and probability 
distributions P(x) show finite size scaling. We note that A/(c) = f™ in (a a ,c) — /™ n (a s ,c). 

accumulate in the same value f x (a x ) where a™ ax = D x and f x = ~{t x — 1)D X . However, for 
probability distribution showing the multifractal scaling there is no accumulation point and 
fx(&x) points shift progressively down as the parameter q is increasing and the parameter q 
approaches D x . This fact is utilized as a simple criterion to recognize which probability dis- 



tributions show either multifractal scaling or FSS [12j. The equality f™ m (a a , c) = f™ m (a s , c) 
is considered to be an attribute that probability distributions P(x) show FSS. To test this 
equality the differences A/(c) denned as A/(c) = /™ n (a a , c) — f™ in (oc s , c) were determined. 
The equality A/(c) = is considered for true if |A/(c)| < 0.10 which reflects numerical 
errors. The differences A/(c) are shown in Fig. H] where the hatched area limits the region 
where the equality is true and thus the avalanche probability distributions P(x) show FSS 
behavior. It is clearly evident that only one value of A/(c) at the density c = 0, is outside 
;he region |A/(0)| > 0.10, and it corresponds to multifractal scaling of the BTW model 
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121 ] . We have no data from the interval of densities < c < 0.01 and thus we may only 



expect that a crossover from multifractal to FSS takes place in this interval. 

The f x (a x ,c) spectra enable us to determine the capacity fractal dimensions D x (c) as 
Dx{c) = a™ ax (c). The results D x (c) for densities < c < 1 are shown in the Fig. [51 
For the BTW model D s (0) = 2.88 ± 0.025 and D a (0) = 2.02 ± 0.025, and for the Manna 
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model D s (l) = 2.77 ± 0.025 and D a (l) = 2.03 ± 0.025 were found. The avalanche area 
capacity fractal dimensions D a (c) are almost constant D a (c) = 2, for any density c, and 
D a {0) = D a (l). In the interval of densities 0.01 < c < 0.15 the avalanche size dimension 
D s (c) is decreasing from D s (0.01) = 2.90 to the value D s (0.15) = 2.78 and is then almost 
constant for c > 0.15, finally D s (0) > D s (l). 

The moment analysis method 12] was used to clarify interesting properties of the scaling 
exponents t x l-+<x{ c ) which are shown in Fig. [TJ The values of the functions f x nin (c) and 
D x (c) ( Fig. E]) are determined from the f x (a x , c) plots. For specific densities c = (the BTW 
model) and c = 1 (the Manna model) /™ n (0) = -0.43 ± 0.05 and f™ in (l) = -0.784 ± 0.05 
were found. Then the scaling exponents are given t x (c) = 1 — f x nm (c)/D x (c) and are shown 
in the Fig. EJ For the density c = 0, it was found T a (0) = 1.213 ± 0.0125. For the densities 
0.01 < c < 0.15, the exponents decrease from r a (0.01) = 1.441 ± 0.0125 and r s (0.01) = 
1.329 ± 0.0125 to the values r a (0.15) = 1.394 ± 0.0125 and r fl (0.15) = 1.299 ± 0.0125, which 
are subsequently constant for c > 0.15. For the density c = 1, they are T a (l) = 1.386±0.0125 
and t s (1) = 1.297 ± 0.0125. These results are similar to those determined directly from the 
distribution functions P(x) ~ x~ Tx (Fig. [T]). 
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FIG. 5: (Color online) The capacity fractal dimensions D x=ajS (c) as functions of the density c. The 
error bars are ±0.025. 
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FIG. 6: (Color online) The scaling exponents were determined using t x (c) = 1 — f x nm (c)/D x (c), 
for the moment analysis all avalanches where a > 50 were taken into account. 

IV. DISCUSSION 

The plots of t X) l vs. 1/ In L and an approximation given by Ecj^ (171) were used to 



extrapolate scaling exponents t Xi £_kx> p, ll7|. Liibeck and Usadel [6|] have analyzed an 
influence of an uncertainty in the determination of the exponents t x> l on the precision of the 
extrapolated exponents t Xj l-><x>- Their results show that this method is not very accurate. 
However, this approximation enables us to make a comparison of our results with previous 
ones. The scaling exponents of the BTW model T a ,L^>oo(0) = 1-26 and t S} l^oo(0) = 1-23 (Fig. 
CP are approximately the same as those found in Ref. [6( (r a = 1.258 and t s = 1.247) using 
the same method. The exponents of the Manna model ^^^(l) = 1.36, and ^^^(l) = 
1.27 are comparable with the previous results, r s L=1024 = 1.28 ±0.02 [2] and with r a = 1.373 
and t s = 1.275, which were found by direct determination of exponents 6] or calculated 
from the moment analysis r a = 1.36 and t s = 1.28 [18|]. The results obtained by the moment 



analysis [12j, / a min (0) = -0.43 ± 0.05 and /™ n (l) = -0.784 ± 0.05, agree well with the 



previous results, a a = -0.391 ± 0.011 and a s = -0.7900 ±0.002 18j. We may conclude that 
the experimental data for two known densities, c = and 1, and data analysis methods 
give approximately the same exponents as were found in previous numerical experiments 



3,0 



18|. 



The scaling exponents defined by Eq. ([6]) [3( and the conditional exponents ^ xy 



can 
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characterize the sandpile models. The theory predicts r s = 1.253 J4] and a few numerical 
experiments show D s ~ 2.7 and D a ~ 2 {4, 10]. The conditional exponents 7™ determined 
directly from the numerical experiments are 7 sa (0) = 1.06 and 7 sa (l) = 1-23 [7]. 

Let us assume that the BTW and Manna models belong to the same universality class. 
Then the scaling exponents r x (c), D x (c) [Eq. ([6] of the model (Eqs. (E])-©) must be 
independent on the density c, i.e. t x (c) = const, and D x (c) = const. This means that 
knowing only the scaling exponents (t x (c), D x (c)), we could not distinguish how many sites 
are toppling by deterministic or stochastic manner [Eq. flS}]. 

We observed that the capacity fractal dimensions D a (c) is constant for any density c, 
D a {c) = 2. The capacity fractal dimension D s (0) = 2.88 is the same as was found in the 



Ref. [12J, D s = 2.86 (determined from the Fig. 
DJl) = 2.77 is higher than the value D ~ 2.7 



(a) in 121]). Our capacity fractal dimension 



101 ]. however it is closer to the D ~ 2.75 



131 ] . In addition, for densities 0.01 < c < 0.1, the scaling exponents t^l^^c) , t x (c) (Figs. 



[Hand El) and D s (c) (Fig. [5]) depend on the density c. These scaling exponents and capacity 
fractal dimension are not constant. They demonstrate that the assumption about a single 
universality class is wrong and thus confirm the existence of different universality classes. 

The conditional scaling exponents j xy [19J can be determined as J xy (c) = (r y (c) — 
l)/(r x (c) — 1) [18]. Substituting the known scaling exponents t x (c) (Fig. E]), we deter- 
mined 7 SO (0.01) = 1.34 and for the Manna model, 7 so (l) = 1.29. We note that the scaling 
exponent t s (0) does not really exist. 

To determine the exact scaling exponents of the probability distribution functions P{x), 
the experimental data must show a power-law dependence given by Eq. ([6]). However, 
the avalanche area size distributions P(a) do not follow exactly power-law distributions 
for densities < c < 0.1 in the whole range of avalanche area sizes, a typical example 
is shown in the Fig. [21 Chessa et al. [lo| found that the area size distribution P(a) of 
the BTW model (c = 0) is not compatible with the FSS hypothesis in the whole range 
of avalanches. However, for large size of avalanches the FSS form must be approached. 
They assume that the scaling in the BTW model needs sub-dominant corrections of the 
form P(x) = (C\X~ Tl + CiX~ T ' 1 + . . .)F(x/x c ) where Cj are nonuniversal constants and 
that these corrections do not determine universality class. The asymptotic scaling behavior 
is determined by the leading power law. We assume that the deviation from a simple 
power-law for densities < c < 0.1 (Sec. IIIII) could be explained by this correction. We 
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observed that the exponents for large avalanches a are larger than the approximate exponents 
(r a L=10 24 = 1-23 in the Fig. [2]) which cover the whole range. As a consequence, the leading 
exponents T a ^(c) for densities < c < 0.1 are higher than the approximate exponents which 
we found (they are not shown in the Figf5]for < c < 0.09). It is evident that the leading 
scaling exponents r a ^ are different and are not constant (Fig. [T]) as in the case of the BTW 
model or the Manna model and thus the model for these densities belongs to a different 
class than the BTW model or the Manna model. 

Divergences from the expected power-law behaviour of the BTW model and a need of 
sub-dominant correction were observed in another inhomogeneous sandpile model 151 ] . Here 
the avalanche dynamic was disturbed by sites which had the second higher threshold. The 
effect was significant for thresholds Eq > 32 and low concentration of such sites [l5|. 

The multifractal properties (Fig. [3]) of the model given by Eqs. (00)-© for the density 
c = (the BTW model), and FSS for the density c = 1 (the Manna mode) agree well with 
the recent results [3]. In addition, the crossover from multifractal to FSS was observed in 
the Fig. HI Our results can only predict that a critical density is expected to be found in 
the interval of densities < c < 0.01 (Figs. [3] and HI). This interval is five times smaller than 

n 

what was found in Ref. [13J where the results are based on the autocorrelation function of 
the avalanche wave time series j^ol ]. 

We assume that divergences from power-law dependences in inhomogeneous conservative 
models, [3] and Eqs. (JTJ)- (EI) , have a common reason which is connected to the crossover 
from multifractal scaling to FSS jl^ . In both models a disorder is induced by deployment 
of disturbing sites. These disturbing sites either increase the short range coupling during 
relaxations in deterministic model [l5|] or introduce the random toppling [Eq. ([5])]. In these 



models toppling imbalance 13|, ]l4| only for a few such sites can change character of waves 
in the models from coherent to more fragmented waves BBS . 121. 

In this study, the multifractal properties of the BTW model which is initially homoge- 
neous, are destroyed at very low concentrations of such disturbing sites. In the opposite 
case, the Manna model shows the FSS and resistance to disturbance caused by presence of 
BTW sites because all significant exponents from Eq. © are approximately constant in a 
broad range of densities 0.15 < c < 1. One possible explanation for this is that the nature 
of the small perturbation of the model is not the same when we perform changes around 
the densities at c = and c = 1. A small perturbation of the dynamical rules of the BTW 
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m papers [7|| or 
scaling P 



model (c = 0) breaks the toppling symmetry [13j and this may explain why the changes 
in the scaling exponents t x (c) and capacity fractal dimension D s (c) are so unexpected. On 
the other hand, for the Manna model (c = 1), decreasing of the density c cannot influence 
the unbalanced toppling symmetry of the Manna model [3]. For sandpile models which 
show FSS this is an expected result and agrees well with the theory % llOj] , where a small 
modification of toppling rules cannot change the scaling exponents. 

We can clearly identify two universality classes which correspond to the classes proposed 
13| | : (a) nondirected models, for density c = (BTW model, the multifractal 
121]). and they show a precise toppling balance 13J and they are sensitive on 
disturbance of avalanche dynamics, (b) random relaxation models, for densities 0.1 < c < 1 
where FSS of P(x) is verified, they are nondirected only on average (Manna two-state model 
c = 1 0]). In these models breaking of the precise toppling balance [3] is observed, the 
scaling exponents are resistant to disturbance of avalanches. The classification for densities 
< c < 0.1 is not so clear. If we follow the proposed classifications then the model is a 
random relaxation model [?| with broken precise toppling balance [3] and it belongs in the 
same class as the Manna model. On the other hand, the scaling exponents differ from the 
Manna model and they are not universal (r x (c) ^ const., D s (c) ^ const.), and the reasons 
of the sub-dominant approximation of area probability distribution functions lo| can play 
an important role. We assume that a new universality class between the BTW (c = 0, 



multifractal scaling) and the Manna (c > 0.5, FSS) classes llJ] could be identified for 



densities < c < 0.1. However, a more detailed study is necessary to verify this classification. 



Our additional arguments to the previous results 



aaaaa 



13 



151 ] show that small 



modifications of the dynamical rules of the model can lead to different universality classes 



what is considered to be unusual from a theoretical standpoint 



io|. 



V. CONCLUSION 



In these comput 
the Manna model 



er simulations multifractal scaling of the BTW model llj and FSS of 



12J were confirmed. In addition, a crossover from multifractal scaling 



to FSS |13j was observed when avalanche dynamics of the BTW model was disturbed by 
Manna sites which were randomly deployed in the lattice, as their density was increased. 
This crossover takes place for a certain density c in the interval < c < 0.01. This interval 



15 



is five times smaller than what was found recently [13j|. The scaling exponents t x (c) and 
the capacity fractal dimension D s (c) are not constant for all densities c which is necessary 



if the models 



HQ 



belong to the same universality class. These result a gree we 



1 with the 



previous conclusions that multifractal properties of the BTW model [5|, [111, [12], toppling 
wave character [3, [si, [q] and precise toppling balance 13, 14 1 are important properties for 
solving the universality issues. 

An open question remains about how to characterize the universality class for densities 
0.01 < c < 0.1, where the scaling exponents are not universal (t x (c) ^ const, and D s (c)j^ 
const.) and in addition, the avalanche probability distributions P{a) do not show exact 
power-law behavior since the sub-dominant corrections of P{a) lOj are important. In this 

7j and to the 



interval of densities c, our model belongs to the random relaxation models 



models with unbalanced toppling sites [13], Q , however, its scaling exponents are not equal 
to the exponents of the Manna model. 

Based on the previous findings k], Q] and our results we assume that the avalanche 
dynamics of undirected conservative models, in which some of the probability distribution 
functions show a multifractal scaling (the BTW model), is disturbed by suitable toppling 
rules which are different from the two-state Manna model (for example a stochastic four-state 



Manna model [7j, |9|), then a local manner for the energy distribution during the relaxation 
can be important and can change the scaling exponents. However, the models which show 
the FSS for all probability distribution functions (the Manna model) are not sensitive to the 
details of the toppling rules and are consistent with theoretical predictions % llO ] - 
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